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Abstract 

I report large-scale Monte Carlo studies of a one-dimensional height-restricted stochastic sand- 
pile using the quasistationary simulation method. Results for systems of up to 50 000 sites yield 
estimates for critical exponents that differ significantly from those obtained using smaller systems, 
but are consistent with recent predictions derived from a Langevin equation for stochastic sandpiles 
[Ramasco et al., Phys. Rev. E69, 045105(R) (2004)]. This suggests that apparent violations of 
universality in one-dimensional sandpiles are due to strong corrections to scaling and finite-size 
effects. 
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I. INTRODUCTION 

Sandpile models are the prime example of self-organized criticality (SOC) , a control 
mechanism that forces a system with an absorbing-state phase transition to its critical 
point P, leading to scale invariance in the apparent absence of parameters Q]. SOC 
in a slowly-driven sandpile corresponds to an absorbing-state phase transition in a model 
having the same local dynamics, but a fixed number of particles P, 0, Q, 0, 0, 13 • The 
latter class of models is usually designated as fixed-energy sandpiles (FES) or conserved 
sandpiles. Continuous absorbing-state phase transitions characterized by a nonconserved 
order parameter (activity density) coupled to a conserved field that does not diffuse in the 
absence of activity, are expected to define a universality class This class, referred to as 
C-DP (that is, a model-C version, in the sense of Halperin and Hohenberg [12J, of directed 
percolation, or DP), appears to be distinct from that of directed percolation [13] . 

In recent years considerable progress has been made in characterizing the critical proper- 
ties of conserved stochastic sandpiles, although no complete, reliable theory is yet at hand. 
As is often the case in critical phenomena, theoretical understanding of scaling and univer- 
sality rests on the analysis of a continuum field theory or Langevin equation (a nonlinear 
stochastic partial differential equation) that reproduces the phase diagram and captures the 
fundamental symmetries and conservation laws of the system. Important steps in this direc- 
tion are the recent numerical studies of a Langevin equation Q, Q for C-DP. (The latter 
which appears to incorporate the essential aspects of stochastic sandpiles.) The critical ex- 
ponent values reported in Ref. are in good agreement with simulations of conserved 
lattice gas (CLG) models 19, |^, which exhibit the same symmetries and conservation laws 
as stochastic sandpiles. 

The Langevin equation exponents are also consistent with the best available estimates for 
stochastic sandpiles in two dimensions [13j, with the exception of the exponent 9 governing 
the initial decay of the order parameter. (The discrepancy regarding 9 likely reflects strong 
corrections to short-time scaling in sandpiles, due to long memory effects associated with 
initial density fluctuations |l^ .) Pending a better understanding of this question, it appears 
that stochastic sandpiles are consistent with C-DP in two dimensions. In the one-dimensional 
case, however, there is a significant discrepancy between the Langevin equation results and 
those for sandpile models. 
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Specifically, analysis of the Langevin equation for C-DP yields, in one dimension, the 
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order-parameter critical exponent value (3 = 0.28(2), while previous studies 3, 
of stochastic sandpiles furnish values near 0.40 for this exponent. There are also smaller 
discrepancies for other critical exponents. If this discrepancy were to persist, one would 
be forced to conclude that the proposed Langevin equation misses some essential aspect 
of sandpiles (at least in the one-dimensional case), or that not all models with the same 
symmetries and conserved quantities belong to the same universality class. In an effort 
;o clarify the situation, I apply the recently devised quasistationary simulation method 
21, 22, 2^ to the restricted-height sandpile introduced in Ref. flfij ]. 

The balance of this paper is organized as follows. In Sec. II I define the model and 
summarize the simulation method. Numerical results are analyzed in Sec. Ill, and in Sec. 



IV I discuss the findings in the context of universality. 



II. MODEL 



16j|. The system, a 



24j, is defined on 



I study the "independent" version of the model introduced in Ref. 
continuous-time, restricted-height version of Manna's stochastic sandpile 
a ring of L sites. The configuration is specified by the number of particles, Zi = 0, 1, or 2, at 
each site i. Sites with Zi = 2 are active, and have a toppling rate of unity. The continuous- 
time Markovian dynamics consists of a series of toppling events at individual sites. When 
site % topples, two particles attempt to move randomly (and independently) to either i — 1 
or i + 1. (The two particles may both try to jump to the same neighbor.) Each particle 
transfer is accepted so long as it does not lead to a site having more than two particles. (If 
the target site is already doubly occupied the particle does not move. Thus an attempt to 
send two particles from site j to site k, with Zk = 1, results in Zk = 2 and Zj = 1.) The next 
site to topple is chosen at random from a list of active sites, which is updated following each 
event. The time increment associated with each toppling is At = 1/Na, where is the 
number of active sites just prior to the event. 

Any configuration devoid of doubly occupied sites is absorbing. Although absorbing 
configurations exist for particle densities p = N/L < 1, the critical value p c (above which 
activity continues indefinitely) appears to be strictly less than unity. In Ref. Q| the 
model was studied in the site and pair mean-field approximation (which yield a continuous 
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phase transition at p c = 0.5 and 0.75, respectively, in one dimension), and via Monte Carlo 
simulation using system sizes of up to 5000 sites. The latter yield the estimates p c = 
0.92965(3), /3/v± = 0.247(2), z = v\\/v L = 1.45(3) and (3 = 0.412(4). A similar value, 
(3 = 0.42(1), was obtained in Ref. |l7[ using a series of cluster approximations (of up to 11 
sites), combined with Suzuki's coherent anomaly analysis 

The studies reported here employ the quasistationary (QS) simulation method, which, 
due to increased efficiency in the critical region, permits a tenfold increase in the system 
size as compared to Ref. The QS method, described in detail in [21], provides a just 

sampling of asymptotic (long-time) properties, conditioned on survival. In practice this is 
accomplished by maintaining (and gradually updating) a set of configurations visited during 
the evolution; when a transition to the absorbing state is imminent the system is instead 
placed in one of the saved configurations. Otherwise the evolution is exactly that of a 
"standard" simulation algorithm such as used in Ref. 



III. SIMULATION RESULTS 



I performed two sets of studies using the QS method. The first is used to determine the 
QS order parameter (defined as the faction p of active sites), the moment ratio m = {p 2 )/p 2 , 
and the mean lifetime r of the quasistationary state, in the immediate vicinity of the critical 
point pc, for system sizes L = 1000, 2000, 5000, 10 000, 20 000 and 50 000. (The QS lifetime 
is taken as the mean number of time steps between successive attempts to visit the absorbing 
state.) A second set of simulations is used to study the supercritical regime (p > p c ) for 
system sizes L = 10 000, 20 000 and 50 000. (For p substantially larger than p c , the lifetime 
is much larger than the simulation time, so that the system never visits the absorbing state, 
and the QS method becomes identical to a standard simulation.) 

Each realization of the process is run for 10 9 time steps; averages are taken in the QS 
regime, which necessitates discarding an initial transient that ranges from 10 6 time steps 
(for L = 1000) to 10 8 time steps (for L = 50 000). The number of saved configurations 
ranges from 1000 (for L = 1000) to 400 (for L = 50 000). The list updating probability p rep 
ranges from 10 -3 (for L = 1000) to 5 x 10 -6 (for L = 50 000). During the initial relaxation 
period P „ p is increased by a factor of ten to erase the rnemory of the tnhia, configuration. 

I first discuss the studies focusing on the critical region. As in Ref. |16J], I study, for each 
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system size, a series of particle number values N, chosen so that p = N/L lies immediately 
above or below p c . Since the particle density can only be varied in steps of 1/L, estimates 
for properties at intermediate values of p are obtained via interpolation. The results of the 
QS simulations were found to agree, to within uncertainty, with the corresponding results 
of conventional simulations |l6j], for L = 1000, 2000 and 5000. The criterion for criticality is 
power-law dependence of p and r on system size, i.e., the familiar relations p ~ L~~PI V± and 
t ~ L z , and constancy of the moment ratio m with L. The most sensitive indicator turns 
out to be the order parameter p. Using the data for system sizes 5000 - 50 000, I rule out 
p values that yield a statistically significant curvature of the graph of In p versus In L. This 
results in the estimate p c = 0.929780(7). (For the remainder of the analysis p c is fixed at 
this value and is no longer available as an adjustable parameter.) The associated exponent 
is (3/u± = 0.213(6), where the uncertainty represents a contribution (±0.005) due to the 
uncertainty in p c and a small additional uncertainty in the linear fit to the data. Simulation 
results for p as a function of L, for various densities near p c , are shown in Fig. 1; curvature 
of the plots for off-critical values is evident in the inset. 

The data for the QS lifetime r furnish a similar but somewhat less precise estimate, 
p c = 0.929777(17). Fitting the data for L = 5000 - 50 000, using the p c interval obtained 
from the analysis of p, I find z = 1.50(4). The moment ratio m is also useful for setting 
limits on p c . As shown in Fig. 2, this quantity appears to grow with system size for p < p c 
and vice-versa; we may rule out the values 0.92976 and 0.92980 on this basis. The moment 
ratio data yield m c = 1.142(8). The main contribution to the uncertainties in z and m is 
again due to the uncertainty in p c . 

The present estimate for p c is significantly greater than that found in Ref . , although 
the difference amounts to about 0.01%. The results for the exponent z are consistent, but the 
present study yields a substantially (16%) lower estimate for /3/u± than reported previously. 
The present result for m c is also substantially lower than the value 1.1596(4) reported in 
Ref. These differences highlight the strong finite-size corrections affecting stochastic 

sandpiles. 

I turn now to the results for the order parameter in the supercritical regime. Fig. 3 shows 
that the data for system sizes 10 000, 20 000 and 50 000 are well converged for A = p — p c > 
10~ 3 , that is, finite-size effects are only present nearer the critical point. Evidently, the 
data are not consistent with a simple power law of the form p ~ A 13 . Indeed this departure 
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from the familiar behavior of the order parameter was already noted (with data for smaller 
systems) in Ref. In the latter work the power law was "restored" by introducing a 

size-dependent critical density p c (L) ~ p coo — Const/ L 1 ^ u± , leading to a series of estimates 
for the critical exponent (3 that increase systematically with L, apparently converging to 
(3 = 0.412(4). With the present data, which are converged over a broader range of A values, 
I find that shifting the critical value does not lead to an apparent power law. 

One is therefore left to conclude that either the order parameter does not obey power-law 
scaling, or that there are unusually strong corrections to scaling. Including a correction to 
scaling term, one has 



P 



\l + AA 13 ' 



(1) 



so that there are now three adjustable parameters, (3, f3' and A. Even with a reasonably 
large number of data points (18 for L = 10 000), this induces a huge range of variation in 
the exponent (3. Decent fits can be obtained with values as low as (3 = 0.1 and as large as 
0.3. 

To resolve this difficulty I return to the data in the immediate vicinity of p c . These 
data can be used to determine the correlation length exponent v± in the following manner. 
Finite-size scaling implies that for p ~ p c , the moment ratio obeys the relation 



m(A,L) ~ J r rn (L 1 / u± A). 
where T m is a scaling function. This implies that 



(2) 



dm 



dp 
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Moreover, the finite-size expression p = L-M v ±F p (L 1 / v± A) implies that 



dlnp 



dp 



oc 
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and similarly for the derivative of lnr at the critical point. The derivatives are evaluated 
numerically as follows. For each value of L studied, data for five values of p clustered around 
p c are fit with a cubic polynomial; the derivative of the polynomial is then evaluated at p c . 
The resulting derivatives are plotted in Fig. 4; clean power laws are observed, leading to 
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u± = 1.362(7), 1.323(14) and 1.372(21), using the data for lnp, m and lnr, respectively. 
Pooling these results yields the estimate u± = 1.355(18). Then, using the values for (3/v± 
and z reported above, I find (3 = 0.289(12) and u\\ = 2.03(8). 

Using this value for (3, the data for the order parameter in the supercritical regime can be 
fit using the correction to scaling form, Eq. Q with parameters (3' = 0.446 and A = 1.3505. 
For A = 0.1, the correction term AA@' in Eq. ^ is 0.48, showing that there are sizeable 
deviations from a pure power law. It is usual to verify scaling by seeking a data collapse, 
plotting p* = L@/ U± p versus A* = L 1 ^ U± A. For A > 0.001 the order parameter does not 
follow a pure power law and so the data cannot collapse. It is nevertheless of interest to 
construct such a scaling plot (Fig. 5). Although the data do not collapse over most of 
the range, they do collapse in the interval — 1 < A < 1. A linear fit to the data in this 
interval yields a slope of 0.27(1). This is close to the (3 value obtained from the finite-size 
scaling analysis, suggesting that simple scaling is restricted to a narrow interval very near 
the critical point. 



IV. DISCUSSION 



A study of the one- dimensional restricted-height stochastic sandpile using quasistationary 
simulations permits study of systems an order of magnitude larger than previously studied, 
and yields critical properties different than those obtained previously. In the case of the 
critical density, the small change (about 0.01%) from the previous estimate may be attributed 
to finite-size effects, which are known to affect sandpile models strongly. 

Of greater concern are critical exponent values, since they define the universality class of 
the model. Since there is every reason (based on symmetry considerations) to expect the re- 
stricted sandpile to belong to the same universality class as the unrestricted version (indeed, 
this seems well established in two dimensions ^f|), I collect, in Table I, critical exponent val- 
ues from various studies of stochastic sandpiles, C-DP and the conserved threshold transfer 
process (CTTP), also expected to belong to the same class. 

The overall conclusion from Table I is that studies using smaller lattices yield values 
in the range 0.38 - 0.42 for the ex pon ent (3 (Ref. is however an exception), and that 
the large-scale simulation of Ref. |2o|. the numerical study of the C-DP field theory 3| 
and the present work yield a consistent set of results, with (3 ~ 0.29. (A similar value has 
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been found for a modified conserved lattice gas model Although the system size (4000 

sites) used in the field theory simulations is not large, one should note that each 'site' in 
such a simulation may represent a region comprising many lattice sites in the original model. 
Compared with the earlier sandpile simulations, the distinctive feature of the present work 
may not be system size, but the fact that here the exponent (3 is determined via finite-size 
scaling at the critical point, rather than from the usual analysis of the order parameter in 
the supercritical regime. Indeed, it is easy to see from Fig. 3 that data for A = p — p c in 
the range 10 -3 - 10 _1 will yield larger estimates for (3. (The same observation applies to the 
CAM analysis ^1, which essentially probes the shape of the function p(A) at some distance 
from the critical point A = 0.) I observe a simple power-law behavior, and data collapse for 
various lattice sizes, only in a restricted range of the scaling variable A* = L 1 ^ V± A. 

Also included in Table I are exponent values for one- dimensional directed percolation 
. The values obtained in Refs. [3| and [2(| , as well as in the present work, are not very 
different from those of DP. A clear difference from DP scaling was however demonstrated 
in Ref. Q|, where the initial decay exponent for one-dimensional C-DP is found to be 
6 = 0.125(2), as opposed to 0.1595(1) for DP. The rather substantial differences found here 
in P/v±, and in the moment ratio m (1.142(8) for the restricted sandpile compared with 
1.1736(1) for DP 2l|). lend further support to the conclusion that the C-DP/stochastic 
sandpile universality class is distinct from that of directed percolation, as is evidently the 
case in two dimensions. (This despite the result j^J, that when suitably modified to include 
'sticky grains', sandpiles fall generically in the DP class.) 

In summary, I have applied the quasistationary simulation method to a one-dimensional 
restricted-height stochastic sandpile, and obtained results consistent with recent studies of 
C-DP. This supports the assertion that the latter class includes stochastic sandpiles, as 
would be expected on the basis of symmetry and conservation laws. 
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TABLE I: Summary of exponent values for one-dimensional models in the C-DP universality class. 
Lmax denotes the largest system size studied. Abbreviations: CAM: coherent anomaly method; FT: 
field theory. 
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FIGURE CAPTIONS 



FIG. 1. Stationary order parameter versus system size for particle densities (bottom to top) 
p = 0.92977, 0.92978 and 0.92979. Inset: lnL a213 p versus lnL for the same set of particle 
densities. 

FIG. 2. Moment ratio m versus system size for particle densities (top to bottom) p = 0.92976, 
0.92978 and 0.92980. 

FIG. 3. Stationary order parameter versus A = p — p c for system sizes (top to bottom) 
L = 10 4 , 2 x 10 4 and 5 x 10 4 . 

FIG. 4. Derivatives of (lower to upper) lnr, hip and m with respect to particle density, 
evaluated at p C) versus system size. The slope of the straight line is 0.734. 

FIG. 5. Scaled density p* versus scaled distance from critical point A*, as defined in text. 
System sizes: 10 4 (open squares); 2 x 10 4 (filled squares); 5 x 10 4 (diamonds). 
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